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Abstract. The one-parametric Wang-Landau (WL) method is implemented together with an extrapolation 
scheme to yield approximations of the two-dimensional (exchange-energy, field-energy) density of states 
(DOS) of the 3D bimodal random-field Ising model (RFIM). The present approach generalizes our earlier 
WL implementations, by handling the final stage of the WL process as an entropic sampling scheme, 
appropriate for the recording of the required two-parametric histograms. We test the accuracy of the 
proposed extrapolation scheme and then apply it to study the size-shift behavior of the phase diagram 
of the 3D bimodal RFIM. We present a finite-size converging approach and a well-behaved sequence of 
estimates for the critical disorder strength. Their asymptotic shift- behavior yields the critical disorder 
strength and the associated correlation length's exponent, in agreement with previous estimates from 
ground-state studies of the model. 

PACS. PACS. 05.50+q Lattice theory and statistics (Ising, Potts, etc.) - 64.60.Fr Equilibrium properties 
near critical points, critical exponents - 75.10.Nr Spin-glass and other random models 

1 Introduction model is 

n = -j s,s,-hY,h^s,, (1) 

The RFIM |ll2l3l4l5l()l7l8l9ll()lllll2li:^ll4ll5j has been ex- <''^> 

where Si are Ising spins, J > is the nearest- neighbors 

tensively studied both because of its interest as a simple 

ferromagnetic interaction, and hi are independent quenched 

frustrated system and because of its relevance to experi- 

random-fields (RF's) obtained here from a bimodal distri- 

ments |16|17|18|19|20|21j . The Hamiltonian describing the 

bution of the form 

e-mail: nfytas@phys.uoa.gr 

^ e-mail: amalakis@phys.uoa.gr P{hi) = -[S{hi - 1) + S{h., + 1)]. (2) 
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h is the disorder strength, also caUed randomness, of the 
system. Various RF probabihty distributions, such as the 
Gaussian, the wide bimodal distribution (with a Gaussian 
width), and the above bimodal distribution [equation ([2])] 
have been considered |22I23I24I25I26I27I28I29I30I31I32I33| . 
As it is well known, the existence of an ordered ferro- 
magnetic phase for the RFIM, at low-temperature and 
weak-disorder, follows from the seminal discussion of Imry 
and Ma [Ij , when the space dimension is greater than two 
(D > 2). This has provided us with a general qualitative 
agreement on the sketch of the phase boundary separat- 
ing the ordered ferromagnetic (F) phase from the high- 
temperature (strong-disorder) paramagnetic (P) phase. 
The phase boundary separates the two phases of the model 
and intersects the randomness axis at the critical value of 
the disorder strength, denoted hereafter as he- Such qual- 
itative sketching has been commonly used in most pa- 
pers for the RFIM |25|31|34|35|36 j and close form quan- 
titative expressions are also known from the early mean- 
field calculations [37> However, it is generally true that 
the quantitative aspects of phase diagrams produced by 
mean-field treatments are very poor approximations. This 
applies also for the bimodal RFIM, for which, with the 
exception of the estimation of he from ground-state cal- 
culations [28129130] , a reliable approximation of the phase 
diagram is still lacking. Furthermore, despite the 30 years 
of theoretical and experimental study the nature and scal- 
ing features of the transition of the RFIM are not yet well 
understood [38139140] . Nowadays, it is generally believed 
that the transition from the ordered to the disordered 
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phase is continuous, governed by the zero-temperature 
random fixed-point [719111] , but a complete set of values 
of the critical exponents fulfilling scaling relations has not 
been established, despite the fact that several bounds [H] 
and further inequalities [8142] for the critical exponents 
have been proposed, together with modified scaling re- 
lations [13]. It is also now quite clear that, the finite- 
size behavior of the system is obscured by strong and 
complex finite-size effects, involving the violation of self- 
averaging [36I44I45I46I47I48I49I5D] . In particular the issue 
of the order of the transition (first-order or continuous) 
has regained much interest after the recent observations 
of first-order-like features at the strong-disorder regime 
for both the bimodal |5T] and the Gaussian RF distribu- 
tions [52153] . 

This work presents a careful and systematic numeri- 
cal approach to the phase boundary of the bimodal RFIM 
in the low-temperature regime. The numerical approach, 
presented below, is a proposal that may be also useful 
to the study of other systems with complex energy land- 
scapes, such as general random systems, spin glasses, pro- 
teins, and others. From our simulations, corresponding to 
systems with linear sizes L in the range L = 4 — 32, we 
perform a finite-size scaling analysis leading also to a re- 
fined value of the critical disorder strength he, in good 
agreement with the estimates obtained via the above men- 
tioned ground-state techniques. We implement a novel ap- 
proach that is based on the idea of entropic sampling in 
restricted energy spaces [54155] together with a reliable 
extrapolation scheme and we produce accurate numeri- 
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cal data at the strong-disorder regime. Our analysis of the study of the RFIM at the strong-disorder regime. Our 
the low-temperature part of the phase diagram provides scheme will be outlined and tested in this Section for the 
us with a qualitative and also quantitative description of 3D bimodal RFIM and it is hoped that it will provide 
the phase diagram of the model, also at low values of the a convenient and fast simulation tool for studying other 
disorder strength, and produces good estimates for the similar disordered or complex systems. In effect, we shall 
critical disorder strength and the correlation length's ex- use our earlier idea of the entropic implementation of the 
ponent, in very good agreement with those from previous WL algorithm ^55J, to produce a faithful approximation 
zero-temperature studies of the model. of the exchange-field two-parametric DOS of the RFIM in 

The rest of the paper is laid out as follows. In the next an appropriate neighborhood of the disorder strength. 
Section we describe the numerical route implemented. In 

The WL algorithm 'B^ is one of the most refreshing 

Section[3]we present in detail the low-temperature aspects 

improvements in Monte Carlo simulation schemes and has 

of the phase diagram of the model. Finally, we summarize 

been already applied to a broad spectrum of interesting 

our conclusions in Section [H 

problems in statistical mechanics and biophysics [63] . Sev- 
eral implementations of the WL sampling technique have 

2 Numerical Approach 

been carried out by many authors |51l52l53l63l64l65l66l67l68l69l70l71l7lj 

There exist two distinct kinds of purely numerical ap- and the present approach may be also seen also as a fur- 

proaches to the RFIM. The first approach utilizes Monte ther contribution to the growing number of different appli- 

Carlo methods, including predominantly sophisticated sim- cations of the WL method in the study of complex systems 

ulation techniques, such as cluster algorithms and flat- with rough energy landscapes. The original WL method 

histogram approaches, to study finite-temperature prop- has been already applied to the RFIM in previous studies 

erties of the system |22l31l34l43l51l56l57l58l59l60j , while concerning the properties of the system at specified values 

the second approach utilizes graph theoretical algorithms of the disorder strength. Such recent investigations have 

to determine the ground-states and estimate the zero- been presented for the bimodal [51] and also for the Gaus- 
temperature behavior of the RFIM |13l27l28l29l30l32l33l52l53iai| RFIM |52I53| . respectively. The present approach fol- 

This second approach, is grounded on the belief that the lows the implementation of the WL random walk used 

critical behavior of the model is governed by the non triv- already in our recent studies of the RFIM [36148149] . In 

ial RF fixed-point at zero-temperature |7I9I11| . these studies we have carried out the WL random walk 

In this work, we follow a novel numerical approach in a restrictive and more efficient fashion. This restric- 

by combining current advances in simulation techniques, five version, utilizes the so called critical minimum energy 

The proposed approach is well adapted and efficient for subspace (CrMES) technique [54l55j to locate and study 
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finite-size anomalies of systems by carrying out the ran- 
dom walk only in tlie dominant energy subspaces. Gen- 
erally, our finite-size scaling studies have shown that this 
restrictive practice can be followed in systems undergoing 
second-order |54l55l69l70l71j and also first-order transi- 
tions |72I73| . Details and tests of this approach for the 3D 
bimodal RFIM can been found in reference [48] , where the 
thermal properties of the system at the disorder strength 
value h = 2 were studied. 

In a subsequent paper |49j the magnetic properties of 
the RFIM were also considered by using the same restric- 
tive scheme as an entropic sampling method. This simpli- 
fication was introduced and tested for the first time in our 
earlier work [55 on the 2D and 3D Ising models and soon 
after that was used for the investigation and verification of 
some universal properties of the order-parameter distribu- 
tion [69j . According to this we may estimate the magnetic 
properties of the systems by recording the two-parameter 
energy-magnetization {E, M) histograms in the final stage 
(high-levels) of the WL diffusion process. At the end of 
the process the final accurate WL (one-parametric) DOS 
G{E) and the cumulative H{E, M) histograms, are used 
to determine the magnetic properties of the system, by 
forming appropriate microcanonical averages of the order- 
parameter moments |49l55l69l71l72l73j . 

The above description may be seen as a convenient 
way to bypass the requirement of a two-parametric WL 
sampling process and a very similar approach will be im- 
plemented in this paper. We will now be recording, again 
in the high-levels of the WL diffusion process, the cumu- 
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lative (exchange-energy, field-energy) two-parametric his- 
tograms, in order to produce an approximation for the 
two-parametric DOS of the RFIM. At this point, we should 
stress that any multi-parametric WL process is inevitably 
restricted to rather small lattices |62I74I75I76I77| . In fact 
the applications of such multi-parametric methods are sub- 
stantially limited, since besides the immense time and ex- 
cessive memory requirements, they very often face severe 
ergodic and/or convergence problems, depending on both 
the physical system and the algorithmic implementation. 
However, notable examples of such two-parametric stud- 
ies, mainly on 2D systems, discussing also some of the 
above problems, have been carried out in the last 10 years. 
The most recent two-parametric investigation performed 
by Tsai et al. 77J concerns the critical endpoint of the 
2D asymmetric Ising model with two and three-body in- 
teractions on the triangular lattice. This last study re- 
quired several days of computer time and a quite large 
computer memory for the larger lattice size studied, con- 
sisting of = 42 X 42 lattice points. To our knowledge, 
this is also the largest system that has been reported by 
the two-parametric WL algorithm. Certainly, a similar 
two-parametric study is possible, although lacking, for the 
RFIM. However, the correspondingly large 3D system will 
have linear sizes of the order of L — 12, and this will be 
very small for our purposes. It will be seen in the next 
Section, that such lattice sizes are rather small for an ac- 
curate estimation of he of the bimodal RFIM. 

We now proceed to give the details of the present en- 
tropic implementation of the WL approach. Carrying out 
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the WL process at a particular value h of the disorder summing over all pairs giving the particular value of the 

strength, we attempt to generate good approximations of total energy 
the (exchange-energy, field-energy) two-parametric DOS 

G,,,{E)= G{Ej,Eh). (4) 

for the RFIM in a neighborhood of h. An analogous ap- Ej+w Eh=E 

proach was undertaken several years ago, before the in- ^et us further assume an entropic Markov process in which 

vent of the WL method, by Deserno [7S], who used flat- "P™ states are selected from the phase space with prob- 

histogram techniques and also a restricted energy sam- ability Wh{x) depending on the DOS Gh{E), where E is 

phng to locate and study some properties of the tricritical t^t^^ energy of the spin state at the value h of the 



point of the Blume-Capel model [79] on a simple cubic lat- 
tice. The extrapolation scheme described below subjects 



disorder, 

Wh{x) cc [Gh{E)]-\ (5) 



to the following WL process: depending on the lattice size, Then, an approximation of the two-parametric (exchange- 
we use a total of at least = 20 WL iterations, produc- energy, field-energy) DOS of the RFIM in a neighbor- 
ing at each iteration level well-saturated energy-histogram hood of h is provided by the expectation of the observable 
fluctuations 80J and obeying at least the 5% flatness crite- 6ej-h jSEh-Hh 
rion |54I55| . The reduction of the WL modification factor , 

follows the usual rule: fj+i = /i = e [54155162) . Mwh(x) ^^^^j^^ 

H^^XEi Eh) 

and the range jwL > 16 of the WL process is used for ~ Gh{E) TrT^ ^ (6) 



the recording and accumulation of the appropriate energy 
histograms (see definitions below). 

To introduce our notation, let us now conveniently sep- 
arate the Hamiltonian of equation ([1]) of the RFIM as fol- 
lows 

n{x) = -JHjix) - hUhix) = -n.jix) - hHhix), (3) 



where the last equality follows from equation ([?]), using the 
above approximate two-dimensional DOS in place of the 
exact and observing that H'^^\E) = Y.Ej+hEu=E H^^'\Ej,Eh) 
and the double histogram H^'^-^Ej^Eh) is the above sum 
of the observable Sej-.-Hj^e^-Mh- The superscript [h) in 
the quantities in the above equation is only a reminder of 



where x denotes a spin state in phase space and we have the fact that the simulation is performed at the disorder 

set J = 1 , since the behavior of the model depends only on strength value h. It should be noted that this notation 

the ration h/ J . Assuming that the two-dimensional DOS does not mean an /i-dependence, but rather a statisti- 

G{Ej,Efi) in the exchange and field variables i?j ~ Ti..j{x) cal indirect influence on the reliability of the histogram 

and Eh = Ti.h{x) is known, the DOS with respect to the recordings and accordingly on the two-dimensional DOS. 

total energy E = Ti{x) = —Ej — h' Eh corresponding to In our approach the ratio of histograms in the above equa- 

any value h! of the disorder strength, can be deduced by tion by the assumed Markov process, is replaced by 
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the ratio developed during the final high-levels {jwL > 16) 
of the WL process. Denoting these latter histograms by 
H'^liEj, Eh) and h'^\{E) and by Gh{E) the WL DOS, 
as modified at the final level of the process, our final ap- 
proximation reads 

G^,iEj,E,)^Gu{E)^^^i§^. (7) 

The above approximation provides in conjunction with 
the skew summing procedure of equation (jj]) a suitable ex- 
trapolation scheme which can be used for the study of the 
RFIM in the neighborhood of the disorder value /i, used 
for the WL simulation. This extrapolation program will be 
applied in the next Section for the study of the finite-size 
development of the phase diagram of the bimodal RFIM 
at the strong-disorder regime. From our previous studies 
it has been verified |55|69j that the detailed balance con- 
dition is quite well satisfied at the high-levels of the WL 
process and the recording of appropriate histograms pro- 
duces faithful and good approximations. Therefore, it is 
hoped that the proposed extrapolation program will not 
produce systematic errors, besides the expected statistical 
fluctuations. However, for safety reasons, we shall use rel- 
atively small values for the extrapolation shift parameter 
\h — h'\, at most of the order of 7% of the disorder strength 
value, and a rather loose restriction of the energy space in 
which the main WL simulation is performed. In particu- 
lar, in most of our simulations performed at h = 2.25 the 
energy spectrum for the simulation was restricted only 
from the high-energy side, while the entire low-energy 
part of the spectrum down to the ground-state was in- 
cluded (see also the discussion below). For the restriction 
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Fig. 1. Illustration of the effect of the violation of the detailed 
balance in the early WL iteration levels. Details of the shown 
approximations are also given in the text. 

of the high-energy side we used our data from our previ- 
ous study of the model at the value /i = 2. In this way the 
WL sampling extends to all energy values with a signifi- 
cant contribution to the finite-size anomalies of the model 
for all values h > 2. For moderately large lattice sizes 
{L > 12), this practice is not the optimum choice. This 
is because, besides the energy states contributing to the 
range h = 2.1 — 2.4, used in our extrapolation program, 
one simulates also the part of the energy spectrum in the 
neighborhood of the ground-states in which the conver- 
gence of the algorithm is very slow. Thus, for the larger 
lattice sizes, one may avoid this ground-state neighbor- 
hood, as we have done for the sizes L = 26 and L = 32. 

Before moving to the presentation of our results, let 
us end this Section by presenting some tests on the reli- 
ability of the proposed approach. Figure [1] illustrates the 
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accuracy of our practice of using the high- levels of the WL 
process as an entropic sampling method. The curves and 
points shown represent three different approximations of 
the specific heat for a particular RF on a lattice of linear 
size L = 8. The solid line is the directly simulated spe- 
cific heat by the WL method at h — 2.1 and should be 
seen as an almost exact result. The open circle points rep- 
resent an excellent approximation obtained for the value 
h — 2.1 by using a WL simulation aX h — 2 and our 
extrapolation scheme, using the high-WL iteration levels 
{jwL = 16 — 20) for the recording the double (exchange- 
energy, field-energy) histograms. Finally, the dashed line 
shows some quite dramatic distortion effects obtained by 
using the whole {jwL — 1 — 20) WL iteration range for the 
recording of the above two-dimensional energy histograms. 
This is of course an example, showing possible subtle ef- 
fects coming from a significant violation of the detailed 
balance condition in the early WL iteration levels. 

A second test showing now the reliability of our ex- 
trapolation scheme is presented in Figure [2l Here we show 
specific heat curves, in the range h = 2.1 — 2.4, obtained 
by the proposed extrapolation scheme from a WL simu- 
lation performed at h — 2.25, together with the results 
obtained independently via direct WL simulation at the 
corresponding disorder strength values. For values very 
close to h — 2.25, the two different approximations coin- 
cide, and even for the values h = 2.1 and h = 2.4 there 
are only very small deviations, mainly around the peaks. 
The locations of the pseudocritical temperatures are very 
weakly dependent on the extrapolation scheme and are 
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Fig. 2. Specific heat curves for a certain RF of the lattice 
size L = 8. Illustration of the behavior for several values of 
the disorder strength obtained by direct WL simulation (lines) 
and by the extrapolation scheme (points) . 

therefore quite accurately determined by the method. The 
effects on ensemble averages will be expected to be even 
weaker. This is illustrated in our final test concerning the 
pseudocritical temperatures obtained from the ensemble 
average specific heat curve, used in the next Section for 
the description of the phase diagram. The average specific 
heat is defined as usually |58|59j 

1 

where the index q — 1, . . . , Q runs over the number of 
disorder realizations. Figure [3] concludes this Section by a 
comparison of two approximations of the average specific 
heat curve [C]av obtained from an ensemble of Q = 35 
realizations and corresponding to the disorder strength 
value h — 2.2. The solid line is the average curve ob- 
tained by a direct WL simulation at h — 2.2, while the 
dashed line represents the approximation of the extrap- 
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Fig. 3. Average specific heat curve at /i = 2.2, obtained by 
direct WL simulation (solid line) and by extrapolation (dotted 
line), for lattice size L = 8 averaged over Q = 35 RF's. 

olation scheme based on a WL simulation on the same 
ensemble at the value h — 2. Clearly, the locations of the 
two pseudocritical temperatures coincide and the two spe- 
cific heat peaks are in excellent agreement. 

3 Phase Diagram 

We aim here to present a reliable approximation of the 
phase diagram of the 3D bimodal RFIM at the strong- 
disorder regime and provide an accurate estimate for he 
(independent from the ground-state approach). Despite 
the general qualitative agreement between different ap- 
proaches on the phase diagram of the model, the various 
estimations throughout the literature vary in a rather wide 
range. This diversity on the numerical estimation of the 
phase diagram is true for both the Gaussian and the bi- 
modal distributions and is generally reflected in the wide 
range of estimates for he- Thus, for the Gaussian RFIM the 
values for he span the range he — 2.26 — 2.37, despite the 
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fact that these values are mainly estimated via the same 
ground-state approach I13I27I28I30I32I33I34I43I53I6T] . On 

the other hand, there are fewer attempts devoted to the 
estimation of the phase diagram of the bimodal RFIM and 
the corresponding estimates for he, obtained again via the 
ground-state approach, are restricted in a smaller range, 
i.e. he = 2.20 - 2.25 |27I28I30| . Our previous attempt 
to estimate the phase diagram using a high-temperature 
(weak-disorder: h — 0.5 — 2) numerical study yielded an 
overestimation for he, namely he — 2.42(18) How- 
ever, we will show here, that an accurate estimation of 
the phase diagram is possible by a more systematic low- 
temperature (strong-disorder: h > 2) numerical study. In 
this case, we will find a much lower estimate for he that 
agrees favorably with the estimates given above from the 
ground-state approach. Additional good comparisons with 
some phase diagram points, estimated earlier in the litera- 
ture, provide evidence that our final proposal for the phase 
diagram may be a reliable and competent approximation 
for the whole disorder strength range. 

We proceed here to analyze our numerical data at 
the strong-disorder regime. Using our cntropic implemen- 
tation of the WL method and the extrapolation proce- 
dure, outlined in the previous Section, we have gener- 
ated numerical data for the following lattice sizes: L S 
{4, 8, 12, 16, 20, 26, 32}. For lattice sizes in the range L = 
4 — 20 we have simulated 20 RF's, whereas for the larger 
sizes L = 26 and L = 32, 10 realizations of the RF have 
been simulated. For each lattice size and each realiza- 
tion, we performed a WL simulation in an appropriate 
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energy subspace, restricted only from the high-energy end 
and including the entire low-energy spectrum down to the 
ground-state, with the exception of the sizes L = 26 and 
L = 32 for which the very close to the ground-sate energy 
levels were avoided. The WL simulation was performed 
at the disorder strength value h — 2.25 and the accu- 
mulated double (exchange-energy, field-energy) histogram 
was then used to approximate the two-parametric DOS 
[equation ^] and finally, the DOS Gh' {E) [equation dJ)] 
and the thermal properties of the system for various values 
of randomness in a neighborhood of the simulated value 
h — 2.25. In order to construct the average specific heat 
curve (Figure [3]) and to identify via its peak a pseudocrit- 
ical temperature T^.h, representing the ensemble of RF's 
at the particular lattice size, a summation over the real- 
izations was performed, as in equation ([8]). As discussed 
earlier and illustrated in Figures [H - [21 the described ex- 
trapolation scheme provides a reliable approximation of 
the location of the maximum of the average specific heat 
curve. The systematic shift of the individual specific heat 
peaks, shown in Figure [21 for higher values of h, will be 
reflected in the corresponding shifts of the peaks of the 
average specific heat curves, as should be expected, pro- 
viding us the necessary information for the finite-size anal- 
ysis. The locations of all these specific heat peaks, for all 
lattice sizes mentioned above, were calculated from our 
simulation data at h — 2.25, and their extrapolations to 
other neighbor /i-values, for the following set I of disor- 
der values, set I: h' = {2.1,2.15,2.2,2.25,2.3,2.35,2.4}. 
For the lattice size L — 12, an additional entropic WL 
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Fig. 4. Approximations of the phase diagram of the 3D bi- 
modal RFIM. Two fitting attempts are shown. The solid line 
corresponds to the elliptical ansatz (|10|l giving he — 2.215(35), 
while the dashed line to the power-law ansatz (|lip giving 
he = 2.277(49). The range of ground-state estimates for he and 
the zero-field's critical temperature Tc-o ~ 4.51153 are marked 
on the axis. The inset shows the shifting of the pseudocritical 
temperature Tl-h for three values of the disorder strength, i.e. 
h = 2.1, 2.15, and 2.2. 

sampling was carried out at h = 2, using now a larger 
ensemble of Q = 250 RF's. Again, using the extrapolation 
procedure of equations ([7|) and (|3|) the specific heat peaks 
corresponding to the following set II of disorder values 
were located, set II: h' = {1.7, 1.8, 1.9, 2, 2.1, 2.2}. 

Let us attempt now a finite-size analysis using the size- 
shifts of the pseudocritical temperatures of the averaged 
specific heat curves for some particular value of the dis- 
order. The inset of Figure [H illustrates fitting attempts 
of these size-shifts for three values of the disorder. The 
range i = 8 — 32 is used in these fits by assuming the 
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usual power law: 

TL;h = T,,h + ahL-^''"^. (9) 

The critical temperatures Tcji, resulting as limiting values 
of the corresponding pseudocritical temperatures, for the 
attempts shown in the inset of Figure HI are 1.297(237), 
0.894(264), and 0.659(299), for the disorder strengths h = 
2.1, 2.15, and 2.2 respectively. We have excluded, from 
our fitting attempts here, the lattice size L = 4 in order 
to eliminate the influence from the very small L-behavior 
and this practice will be followed and further discussed in 
the sequel. Following the same fitting procedure, again in 
the range L — % — 32, for h — 2.25 we find that the corre- 
sponding critical temperature becomes now negative, i.e. 
?c;2.25 — —0.18. This fact shows that, within our fitting 
scheme, the value of the disorder strength h ~ 2.25 is 
an upper bound for the critical disorder strength. Note- 
worthy, that if we use the range L = 4 — 32 instead, the 
negative sign for the critical temperature will appear at 
the value h — 2.35, which however appears to be a rather 
overestimating bound for the critical randomness. Thus, 
only the three points h — 2.1, h = 2.15, and h = 2.2 
(filled circles) resulting from the fits shown in the inset 
of Figure |4] can be used to approximate the phase dia- 
gram. In order to find one more point of the phase dia- 
gram we shall now also use our earlier numerical data [48] 
(from rather large Q — 500 — 1000 ensembles of RF's) 
for the disorder strength h — 2. Using the above fitting 
practice in the range L = 8 — 32 we find from the gen- 
eral pseudocritical temperature shift behavior the limit- 
ing value Tc-2 = 1.848(188) (open triangle), which is just 
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inside the estimate bounds given in our previous paper 
{Tc-2 = 2.03(18)) using sizes in the range L = 4 - 32 [48] . 

The above four approximate phase diagram points, 
corresponding to the disorder strength values h — 2, 2.1, 
2.15, and 2.2, will be now used to find a phcnomenolog- 
ical representation of the phase diagram of the bimodal 
RFIM. Let us first attempt an elliptical fit using the fol- 
lowing ansatz 




The rescaling temperature factor t in equation ()10|) will 
be handled either as a free-parameter during the fit, or as 
a fixed-parameter using the best known estimate for the 
critical temperature of the zero-field Ising model, namely 
Tc:0 — 4.51153 [81 . The resulting phase diagrams almost 
coincide (see Figure |4] where for clarity reasons only the 
latter case is shown) and are described respectively by the 
following {he, T, x; x^) parameter values, including the 
value of the x^-tcst: (2.212(29), 4.50394(778), 1.862(87); 
- 10--*) and (2.215(35), 4.51153, 1.847(92); - lO""*), re- 
spectively. Thus, our fitting attempts with equation (jlOp 
produce a value for the critical disorder which is very close 
to the estimates obtained from the zero-temperature stud- 
ies of the model )27I28I30] . Furthermore, the fitting using 
the temperature rescaling factor t in equation (|10p as a 
free-parameter produces a fairly good estimate for the crit- 
ical temperature of the zero-field Ising model [5T] . 

As an alternative to the above elliptical fit, we have 
also considered for comparison the following power-law 



N.G. Fytas and A. Malakis: Phase Diagram 

ansatz [36] 

h=hc{^^^^^y. (11) 

The attempt to fit the same data to this law is illustrated 
also in Figure [4] by the dashed line. In this case we find 
a noticeable overestimation of namely he — 2.277(49) 
and a much larger (by a factor of 70) value of of the 
fit. Therefore, we conclude that the elliptical law of equa- 
tion pU|) provides a better representation of the phase 
diagram of the RFIM. Of course, our attempt above aims 
only at a numerical approximation for the main part of the 
diagram and not at the correct asymptotic behavior at its 
ends. For instance, the behavior of the phase diagram at a 
very small neighborhood around the critical temperature 
of the pure system, is expected to be determined by the 
susceptibility exponent 7 of the pure system [37139182] , as 
follows from the phenomenological renormalization argu- 
ments of reference [31] . Accordingly, the slope of the phase 
diagram at this end is expected to behave as 5h ^ {5Ty 
(where 7 = 1.2358 for the pure 3D Ising model [SSJ) and 
not with the exponent 1/2 of the ansatz (tTUj). It appears 
that similar elliptical laws have been also used previously 
by other authors for the Gaussian RFIM [31134] , although 
these were not stated explicitly. 

Finally, we would like to note that we have included in 
Figure H] some more data points for smaller values of the 
disorder strength from previous numerical works. These 
are the data for h — 0.5, 1, and 1.5 (shown by stars in the 
figure) from our previous investigation of the phase dia- 
gram of the model [36| and three more points (open cir- 
cles) estimated by Rieger and Young [58j . These points are 
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close enough to our approximate phase diagram and the 
small deviation comes possibly from the fact that these 
have been estimated, in both cases, by applying equa- 
tion ^ to rather small sizes: L < 2A and L < 16, re- 
spectively (see also the discussion below). 

At this point, let us comment on the significance of our 
notation concerning the shift exponent i>h in equation @ . 
As mentioned earlier, we have tried to avoid the influence 
of the very small L-behavior in our estimates, thus exclud- 
ing from our fitting attempts the data for i = 4. This is 
a compromise followed because in our study (and in ef- 
fect in all finite-temperature studies) a rather restricted 
i-range is available for performing finite-size scaling anal- 
ysis. However, it has been pointed out in reference 36J that 
the estimates based on such restricted ranges should not 
be completely trusted and this may be particularly true for 
the shift exponent Vh. For instance, the range L = 4 — 24 
will produce quite different estimates, for both Tc-h and 
i>h, from those obtained above from the range i = 8 — 32 
and this fact, together with the use of the power-law in 
equation pip , are the two reasons behind our earlier over- 
estimation of he {he — 2.42(18) in reference [35]). The 
general asymptotic behavior of the RFIM follows differ- 
ent complex routes that appear to strongly depend on the 
value of the disorder strength h and different ranges of lat- 
tice sizes may be needed in order to approach the asymp- 
totic behavior for different values of disorder strengths. 
Even the observation of an appreciable disorder strength 
dependence on v^, should be reluctantly identified as a 
possible violation of universality along the phase bound- 
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Fig. 5. Finite-size elliptical phase diagrams for L — 12 us- 
ing two different extrapolation sets of disorder strength values 
(set I (filled triangles) and set II (open circles)) and different 
realizations ensembles. The solid and dashed lines are ellipti- 
cal fits of the form (|12[1 with comparable values of of the 
order of 10~^ giving for the pseudocritical disorder strength 
the values hc-i2 = 2.56(3) and /ic;i2 = 2.56(2), respectively. 
The application of the finite-size version of the power-law 
shown by the dotted line, has a larger value of ~ 10~^ 
and produces an overestimation for the pseudocritical disorder 
strength: /ic;i2 = 2.77(5). 



ary, although this violation of universality is one of the 
strongly supported scenarios in the literature [53] . The vi- 
olation of universality for the case of the 3D RFIM has 
been discussed a few years ago by Sourlas |29j. Equiva- 
lent studies of universality violations have been reported 
also in other glassy systems [53] , reenforcing the view that 
the concept of universality in complex systems is not fully 
clarified. 



We proceed now with an alternative estimation of the 
critical disorder strength. Firstly, let us point out that 
for each value of L, our data can be used to produce a 
finite-size phase diagram. Provided that the phase dia- 
gram points do not decline appreciably from the above 
elliptical law, we may attempt to construct a finite-size 
sequence of diagrams by using the finite-size version of 
equation ([TI 



(12) 



where now the rescaling temperature factor tl may be 
either handled as a free-parameter during the fit or as 
a fixed-parameter at the corresponding zero-field's Ising 
model pseudocritical temperatures taken from Table IV 
of reference [54] . Using this latter choice for tl , Figure [5] 
provides a test of this approach producing two very simi- 
lar phase diagrams for the size L — \2. The two diagrams 
are obtained using the two different sets of phase diagram 
points corresponding to set I and set II of the disorder 
strength values. The first set of points (filled triangles) is 
determined over an ensemble of Q = 20 realizations of the 
RF and corresponds to set I, i.e. simulation at ft. = 2.25 
and suitable extrapolation in the range h — 2.1 — 2.4. 
The other set of points (open circles) is determined over 
a larger ensemble of Q = 250 realizations of the RF and 
corresponds to set II, i.e simulation at ft = 2 and extrapo- 
lation in the range ft — 1.7 — 2.2. The application of the el- 
liptical law (I12|) gives the two very similar phase diagrams 
shown in Figure[5]by the solid and dashed lines for the two 
set of points, respectively. These two diagrams, with com- 
parable values for of the order of 10"^, come together 
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to the same point at T = 0, giving the value ft.c;i2 = 2.56 
for the finite-size {L — 12) critical disorder strength. For 
illustration reasons, we have also included in this figure 
the attempt using the corresponding finite-size version of 
equation for the set of points obtained from the sim- 
ulations at the value h = 2.25 (dotted line). Again the 
quality of the fit is worst for the power-law (x^ = 10^^) 
and produces a clear overestimation for the pseudocriti- 
cal disorder strength of the order of hc-i2 — 2.77(5). The 
comparison between the two finite-size elliptical phase di- 
agrams, corresponding to the two sets of points {h — 2.25 
and h — 2), is on the other hand very convincing. Thus, 
Figure [5] provides strong evidence in favor of our choice 
of using in our simulations for all lattice sizes the strong- 
disorder regime corresponding to the value h — 2.25. In 
particular it shows that our data based on only the Q = 20 
RF's are sufficient for our proposes of estimating the phase 
diagram. 

Figure [H] presents the finite-size elliptical phase dia- 
grams for lattice sizes in the range L = 8 — 32, using 
set I of the disorder strength values. For the lattice size 
L = 4 we have not drawn a finite-size phase diagram, 
since it is quite evident from the corresponding open star 
points in Figure [5] that they decline very early, at about 
the value h = 2.2, from the elliptical law. No such de- 
viation is observed for the other lattice sizes, within the 
set I of disorder values, and this fortifies our choice to 
use the particular set I for these lattice sizes. Of course, 
an attempt to push our approach to even larger lattices 
may require a WL simulation at h = 2.2 and a corre- 
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Fig. 6. Finite-size elliptical phase diagrams for lattice sizes in 
the range L = 8 — 32, using set I of the disorder strength val- 
ues. The drawn lines represent the finite-size elliptical fittings 
according to equation (|12p . in which the rescaling tempera- 
ture factor TL was fixed at the corresponding zero-field's Ising 
model pseudocritical temperatures. 

spending set of somewhat smaller disorder values. The 
drawn lines in Figure [S] represent the finite-size elliptical 
fittings according to equation p2|) . in which the rescal- 
ing temperature factor was fixed at the corresponding 
zero-field's Ising model pseudocritical temperatures. For 
clarity the diagrams using tl as a free-parameter are not 
shown. However, the main frame and the inset of Figure [7] 
illustrate the smoothness of the both fitting schemes and 
reveal a convincing and regular shift-behavior of the finite- 
size critical disorder strengths /ic:L- This behavior allows 
now a finite-size analysis for the estimation of he- The solid 
and dashed lines show good quality fits to the following 
usual shift power-law 

hc.L^hc + bL-^/". (13) 
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Fig. 7. Shift behavior of the finite-size critical disorder 
strengths hc-L- The inset shows the osciUations in the values of 
tl, when this parameter is handled as a free-parameter. Their 
behavior follows the correct trend, approaching the zero-field's 
Tc;o = 4.51153 tSlj (dotted line). 

Thus, the fitting attempts in Figure [7] produce estimates 
for the asymptotic value of the critical disorder strength 
hc^ and the corresponding shift-exponent v. The fitting 
scheme based on the estimates of the pseudocritical dis- 
order strengths h^j, (open circles in Figure [7|), produced 
by fixing the rescahng temperature factor tl at the corre- 
sponding zero-field's Ising model pseudocritical temper- 
atures, i.e. tl = Tl:o, gives he = 2.219(83) and ly — 
1.806(390). Finally, the fitting attempt based on the cor- 
responding hc:L estimates (filled triangles in Figure [7]), 
produced by using tl as a free-parameter, results in a al- 
most identical estimate for the critical disorder, i.e. he = 
2.219(65) but a slightly lower estimate for the shift-exponent 
1/ — 1.640(423). From the inset of Figure [7] we may note 
some oscillations in the values of tl , when this is handled 
as a free-parameter, which however appear to follow the 
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correct trend so that will approach Tc-q (dotted line) 
with increasing lattice size. In both cases, the estimates 
for he compare very well with those obtained above from 
fitting equation (fTU|) in Figure [3] and also with those of 
the ground-state approach |27I28I30] . Despite the devia- 
tion of the two estimates for the shift-exponent and the 
relatively very large variation of ly in the literature (for 
both the Gaussian and bimodal cases), it is of interest to 
compare here the estimate of the second case {ly = 1.64) 
with the estimates 1.67(11) and 1.66(8) of references [25] 
and [30] respectively, obtained by zero-temperature simu- 
lations. 

The above observations provide concrete evidence in 
favor of our present approach. It appears that, this method 
may be capable to produce, if further pushed to larger lat- 
tices, even more accurate estimates for both the critical 
disorder strength and also the T = correlation length 
exponent, assuming that its behavior follows the observed 
shift-behavior of our finite-size projections /ic;L- It is well 
known from the general scaling theory that, even for sim- 
ple models, the equality between the correlation length's 
exponent and the shift exponent is not a necessary con- 
sequence of scaling [85]. Of course, it is a general prac- 
tice to assume that the correlation length behavior can 
be deduced by the shift behavior of appropriate thermo- 
dynamic functions. In our view, the recent strong version 
of the zero-temperature fixed-point scenario by Wu and 
Machta |52I53| . supports the above assumption that the 
finite-size projections /ic:L are appropriate shifting param- 
eters. The thermal states of Wu and Machta (see Figure 
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4 of reference [52]) at temperatures close to the finite- 
size anomalies are strongly correlated to the ground-states 
at disorder strength values close to the zero-temperature 
critical point and this strong correlation may be seen as a 
phenomenological justification of our assumption. 

4 Conclusions 

A numerical approach combining well-known techniques 
has been proposed as a convenient alternative for the study 
of disordered systems. Within this approach, the well- 
known WL algorithm is used, at its final stage, as an cn- 
tropic sampling method, and multi-parametric histograms, 
appropriate for the study of the system, are produced. The 
main advantage of this scheme is that the requirement of 
multi-parametric WL sampling is surpassed and by using 
the DOS, obtained via the WL method, and the accumu- 
lated histogram information, the thermal properties of the 
disordered system may be obtained in a neighborhood of 
the simulated disorder strength value. The numerical tech- 
niques presented in this paper may find further applica- 
tions in the study of critical properties of other challenging 
disordered systems. Via the above approach, we have stud- 
ied the general size-shift behavior of the low-temperature 
part of the phase diagram of the 3D bimodal RFIM. Our 
detailed analysis provided an overall reliable representa- 
tion of the main part of the phase diagram, yielding ac- 
curate estimates for the critical disorder strength. These 
estimates are in agreement with those from previous zero- 
temperature studies of the model including the estimates 
for the correlation length's exponent. 
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As a closing remark, we would like to mention that, 
using our WL DOS's - for some typical RF realizations, 
at the simulated disorder strength value h = 2.25 - we 
have also observed, for the larger sizes studied, first-order- 
like double peaks in the energy probability densities, in 
agreement with the recent observations of Hernandez and 
Ceva [5TI, and Wu and Machta [52l53j . mentioned in the 
introduction. This main issue appears to be still a mat- 
ter of controversy and we are currently carrying out fur- 
ther research in order to clarify the persistence (or not) 
of such first-order-like characteristics in the asymptotic 
limit. However, the full resolution of this aspect requires 
an understanding of the complex finite-size effects of the 
RFIM at the strong-disorder regime and substantial com- 
puter resources to be devoted for the simulation of large 
ensembles of RF realizations in a convenient neighborhood 
of disorder strength values. 
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